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Abstract 

Integrating biological information from different sources to understand cellular pro- 
cesses is an important problem in systems biology. We use data from mRNA expression 
arrays and chemical kinetics to formulate a metabolic model relevant to K562 ery- 
throleukemia cells. MAP kinase pathway activation alters the expression of metabolic 
enzymes in K562 cells. Our array data show changes in expression of lactate dehydro- 
genase (LDH) isoforms after treatment with phorbol 12-myristate 13-acetate (PMA), 
which activates MAP kinase signaling. We model the change in lactate production 
which occurs when the MAP kinase pathway is activated, using a non-equilibrium, 
chemical-kinetic model of homolactic fermentation. In particular, we examine the role 
of LDH isoforms, which catalyze the conversion of pyruvate to lactate. Changes in the 
isoform ratio are not the primary determinant of the production of lactate. Rather, 
the total concentration of LDH controls the lactate concentration. 
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Introduction 

Modeling of cellular metabolism has a long history of important contributions to 
biology. Approaches include kinetic modeling, metabolic control analysis, flux bal- 
ance analysis, and metabolic network analysis ( |Stephanopoulos et al, 1998| ). A new 
era of research in metabohsm is now possible, because large-scale expression studies 
can determine levels of many metabolites ( |Goodac~ et al., 2004t|Fan et al., 2004|) and 



metabolic enzymes dFerea et a/.,T999l|Kar7^ al., 1999D . In this paper we use metabolic 



enzyme expression data to guide metabolic modeling, with a focus on small but signifi- 
cant changes in mRNA abundance. Our goal is to understand how biologically realistic 
changes in mRNA abundance of metabolic enzymes affect cellular metabolism. Previ- 
ous work integrating expression data with metabolic modeling has been done in yeast 
( Akesson et al., 2004D , but not, to our knowledge, in mammalian systems. 



We focus on glycolysis, an essential ATP-producing metabolic pathway. The initial 
reactions of glycolysis break down glucose into pyruvate. Pyruvate can feed into ei- 
ther the citric acid cycle (aerobic metabolism) or homolactic fermentation (anaerobic 
metabolism) ( |Voet fc Voet, 2004D . The reactions involving pyruvate therefore control 
this important metabolic branch point. Homolactic fermentation is catalyzed by lactate 
dehydrogenase (LDH) in a compulsory-order, ternary reaction (Borgmann et al., 1975 1 . 



LDH reversibly converts pyruvate and NADH into lactate and NAD-I-. The isozymes 
of LDH are tetramers formed from two types of monomers (a third isoform is usually 
germ-line specific, but can be expressed in cancers ( IKoslowski et al., 2002D ). The two 



isoforms are labeled H (heart) and M (muscle), and their ratio varies between cell 
types. The LDH isoform ratio has been proposed to indicate the metabolic state of 
cells: it is believed that the M isoform favors lactate production while the H isoform 
favors pyruvate production ( |Boyer et al, 1963 Stambaugh fc Post, 1966| |Boyer, 1975| 



Voet fc Voet, 2004| ). In this framework, the LDH isoform ratio can serve as an indicator 



of the relative flux through aerobic/anaerobic gycolytic pathways. 

Here we use a mathematical model of homolactic fermentation to study the con- 
nections between growth-factor signaling and metabolism. It has been known since the 
work of Warburg that carcinogenesis is accompanied by changes in cellular metabolism 
( Stubbs et al, 2003 Griffiths et al, 2002[ [Dang fc Semenza, 1999 1. In particular, tu- 



mors typically favor anaerobic metabolism, resulting in higher lactate production rela- 
tive to non-cancerous cells dWalenta et al., 2004| [Newell et al., 19931 [Warburg, 19561 ). 
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Although inhibition of glycolysis can kill tumor cells ( Munoz-Pinedo et al., 2003 1, the 
connections between carcinogenesis and metabolic alterations are not fully understood 
( Fan et al, 2004 1 . However, intriguing connections between metabolic enzymes and 
cancer have been demonstrated ( Kondoh et al, 2005) Kim et al., 2004[ [Mazurek fc Eigenbrodt, 2003 1 . 
In particular, LDH expression is altered in many tumors dWalenta et al., 2004| |U nwm et al., 20031 
[Maekawa et al, 2003 1 and cancer cell models ( |Li et al., 20041 |K aran et al., 20021 |L ewis et al, 2000| ). 
High tumor LDH levels have been shown to correlate with poor prognosis in lung cancer 
patients ( |Koukourakis et al, 2003 1. 

In this study, we focus on changes in LDH expression induced by the mitogen- 
activated protein (MAP) kinase pathway. The MAP kinase cascade is important in 
cell growth, differentiation, and survival, and alterations of MAP kinase signaling have 
been found in many cancers ( Lewis et al, 1998 ). The signal is transduced by a series 
of phosphorylation reactions: MAP kinase proteins phosphorylate and thereby acti- 
vate their downstream targets. The pathway includes the MAP kinase proteins ERK 
1 and 2 and their upstream activators, the MAP kinase kinases MKK 1 and 2. Re- 
cent work has found connections between MAP kinase signaling and metabolism. For 
example, increased expression of LDH-H in human tumors may occur in part because 
the transcription factor MYC, a downstream target of the MAP kinase pathway, tran- 
scriptionally up-regulates the LDH-H gene ( Jungmann et al, 1998|l^him et al., 1997D . 
Other genes involved in glycolysis are also affected by MYC dUsthus et al, 20001 ). Ac- 
tivation of the MAP kinase pathway has been shown to increase LDH activity, glucose 
uptake, and lactate production dRiera et al., 20031 IPapas et al., 1999D . 

We studied mRNA expression in K562 erythroleukemia cells, a cell line used as 
a model for leukemia. In our experiments, MAP kinase signaling was either [i) ac- 
tivated with phorbol 12-myristate 13-acetate (PMA) or (ii) simultaneously activated 
with PMA and inhibited with U0126, a specific MKK inhibitor. We found small but 
reproducible changes in the expression of LDH isoforms in response to MAP kinase 
pathway activation (figure Q), with no significant changes in other enzymes that cat- 
alyze reactions involving pyruvate. This result suggests that activating the MAP kinase 
pathway alters the relative fiux through aerobic and anaerobic glycolysis in these cells. 
We chose to model the expected changes in cellular lactate prodution to better under- 
stand the connections between signaling and metabolism. We hypothesized that the 
LDH isoform ratio plays an important role in determining cellular lactate levels, as 
suggested previously ( ]Riera et al., 20031 IDang & Semenza, 1999D . 

We formulated a chemical-kinetic model of homolactic fermentation based on in 
vitro biochemistry dBorgmann et al., 1975 ). Our goal was to determine how changes 
in the LDH isoform ratio alter the amount of lactate produced by K562 cells. We 
used the experimentally determined abundance changes as model inputs. The model 
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describes the mass-action kinetics of homolactic fermentation. We included metabolic 
flux terms in the model to describe the connection between homolactic fermentation 
and the larger metabolic network of the cell. The metabolic flux is a constant rate of 
production/consumption of a metabolite through other reactions or transport. Several 
model inputs — the steady-state concentrations of pyruvate, NADH, and NADH — have 
not been measured in K562 cells. Therefore we validated our results with a robustness 
analysis ( |von Dassow et al, 2000| [Barkai fc Leibler, 1997D . 

We present several unexpected flndings. In a preliminary analysis, we examined the 
behavior of each isoform individually. Our results predict that LDH-H produces a larger 
steady-state lactate concentration than an equivalent amount of LDH-M under typical 
cellular conditions. This result is surprising because it disagrees with the statement, 
often found in the literature, that the M isoform favors lactate production while the 
H isoform favors pyruvate production ( |Boyer et al, 1963| Stambaugh &: Post, 1966| 
IBoyer, 19751 |Voet & Voet, 2004D . We discuss the reason for this difference and explain 
why our results are more applicable in vivo. 

Second, we predict a decrease in the steady-state lactate concentration when the 
LDH isoform abundance shifts from control to PMA-treated levels. This result means 
that the H:M isoform ratio alone does not control the lactate concentration. After 
PMA treatment the ratio of LDH-H to LDH-M changes from 1.02 to 1.35 in our ex- 
periments. According to our single-isoform model, an increasing isoform ratio should 
lead to an increase in lactate concentration. This flnding led us to consider separately 
how the isoform ratio and the total abundance of LDH control the lactate concen- 
tration. We demonstrate that while the isoform ratio does affect the production of 
lactate, the experimentally determined total LDH abundance change plays a larger 
role in determining the lactate concentration. 



Methods 



Cell extraction and microarray analysis 

K562 erythroleukemia cells were grown in suspension in 10% FBS/RPMI and 
treated with 10 nM phorbol 12-myristate 13-acetate (PMA) and 20 fiM U0126 (Promega) 
as described previously dSevinsky et al, 2004 1). Cells (7 x 105) were washed twice in 
ice cold phosphate buffered saline, 1 mM EDTA, 1 mM EGTA, and total RNA was 
isolated by TRlzol extraction (Invitrogen). First and second strand cDNA synthesis. 
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in vitro transcription of biotin-labeled cRNA, and fragmentation were carried out fol- 
lowing standard protocols from Affymetrix. Probes were hybridized onto U133 2.0 Plus 
GeneChips (Affymetrix) and processed at the UCHSC Cancer Center Microarray core 
facility. Datasets were corrected for background and normalized using RMA Express 
software. Each condition was analyzed in three independent experiments (figure^. 

Model 

We used a chemical-kinetic model to analyze the effect of isoform switching on the 
non-equilibrium steady state of homolactic fermentation. We assume that the metabo- 
lite and enzyme species are homogeneously distributed in the cytosol. This leads to a 
set of 12 mass-action ordinary differential equations that describe the time evolution of 
metabolite and enzyme concentrations. Four equations govern the metabolites and four 
equations govern each of the LDH isoforms and their related complexes (equations 

Each elementary reaction in the model (figure Ej) follows the law of mass action, 
which results in the following reaction rates 



vi = kiXiCi - k_ie2, 


(la) 


V2 = k2X2e2 - A;_2e3, 


(lb) 


V3 = ^363 - ^-32/1 64, 


(Ic) 


f4 = ^464 - k_4y2ei, 


(Id) 


V5 = k^xifi - A;_5/2, 


(le) 


v& = keX2f2 - k^efs, 


(If) 


vr = kjf^ - k^jyifi, 


(Ig) 


vs = ksfi - k^sy2fi, 


(Ih) 



where xi, X2, yi, and y2 are the concentrations of NAD+, lactate, pyruvate, and 
NADH; ci, 62, 63, and 64 are the concentrations of LDH-H, LDH-H:NAD+, LDH- 
H:NAD-|-:lactate, and LDH-H:NADH; fi, /2, /3, and /4 are the concentrations of LDH- 
M, LDH-M:NAD+, LDH-M:NAD+:lactate, and LDH-M:NADH. The values of the 
kinetic rate constants are shown in figure El 

The equations describing the dynamics of the system can be written compactly in 
terms of the reaction rates. The equations for the metabolites are 



x[ = -vi - V5 - ai, (2a) 

X2 = -V2 - vq - a2, (2b) 

y'l = V3 + VY + Q!3, (2c) 

y2 = Vi + Vs + ^4, (2d) 
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where ai and 02 are the flux of NAD+ and lactate out of the system, and and 04 
are the flux of NADH and pyruvate into the system. 

LDH is an efficient catalyst which is thought to operate near equilibrium under 
many conditions ( |Borgmann et al., 1975D . However, the metabolites in homolactic fer- 
mentation are continually consumed or produced in other reactions or transported into 
and out of the cell. This flux of metabolites means that the system is not in ther- 
modynamic equilibrium. In contrast to an equilibrium system, the steady state of a 
non-equilibrium reaction depends on the reaction mechanism and the total concentra- 
tion of the enzyme that catalyzes it. Here we do not specify a mechanism for this 
metabolic flux but assume that each metabolite is produced or consumed at a constant 
rate. The constant flux of each metabolite in the model is represented by the constant 
terms ai , . . . , 04 in equation (j2)) . 

The heart-isoform complexes obey the equations 



e'l = f4 - vi, 


(3a) 


e'2 = Vi- V2, 


(3b) 


e'^ = V2- V3, 


(3c) 


64 = f3 - Vi, 


(3d) 


obey the equations 




f[=Vs- V5, 


(4a) 




(4b) 


/g = V6 - V7, 


(4c) 


fi = VT- Vs. 


(4d) 



We are primarily interested in the steady-state behavior of the model, which occurs 
only when the fluxes are equal for all four metabolites. In other words, a steady state 
exists if and only if ai = a2 = 0:3 = = a, where a is the constant metabolic flux of 
the model. 



Numerics 



A numerical approach was used to determine the steady-state relationship between 
the metabolite concentrations, the metabolic flux, and the isoform concentrations with 
both isoforms present (flgure EJ. Equations (PI)-© can be integrated numerically, 
allowing the system to approach a steady state from any initial condition. However, 
we specifled the steady-state concentrations of NAD+, NADH, and pyruvate, and 
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determined the corresponding steady-state concentration of lactate (equations 
(Pc|) . and ()2d|) were eliminated from the model). 

The differential equations were integrated using a Runge-Kutta-Fehlberg fifth-order 
method with adaptive stepping ( |Press et al, 1992 ). Newton's method was used to ac- 
celerate convergence to the steady state. If the solution to the constrained differen- 
tial equations was close to a steady state, Newton's method quickly converged to it 
( Press et al, 1992 ). If Newton's method started too far from the steady state it rapidly 
diverged. When this happened, the constrained equations were integrated again (con- 
tinuing from the last solution) for a specified time interval. Integrating for a longer 
time allowed the system to approach closer to the steady state before Newton's method 
was tried again. This hybrid approach worked well for all the conditions we studied. 

We note that Newton's method does not follow the dynamics of the system and 
does not find the correct steady state unless the equations are additionally constrained. 
In particular, it is necessary to explicitly enforce the conservation of the total enzyme 
concentrations (see below). 



Results 



Single-isoform results 

We examined the steady-state production of lactate by a single LDH isoform. Under 
typical cellular conditions our model predicts that LDH-H produces more lactate than 
LDH-M. This result is surprising because many authors state that LDH-M produces lac- 
tate more efficiently than LDH-H ( |Boyer et al., 1963| [Stambaugh fc Post, 19661 jBoyer, 1975| 



Voet fc Voet, 2004D . We expl ain the reason for this difference, which results from dif- 



ferent model assumptions, and argue that our analysis is more experimentally relevant. 



Steady-state lactate production by a single isoform 

Here we are interested in the behavior of the model at steady state. All concentra- 
tions are assumed to be steady-state values unless otherwise stated. Suppose that the 
concentrations of pyruvate, NAD+, and NADH are known. Furthermore, suppose only 
one isoform of LDH is present and its total concentration is known. We can then derive 
an equation that describes the concentration of lactate as a function of the metabolic 
flux a. Recall that when a > lactate and NAD-I- are removed from the system and 
pyruvate and NADH are added to the system. 



7 



The total concentration of either enzyme isoform is constant. This can be shown 
for LDH-H by adding equations (jHajl - pHj) . 

e\ + 62 + + = 0, (5) 

and integrating to obtain 

ei + 62 + 63 + 64 = 6o. (6) 

Here eo is the total concentration of the heart isoform. This relation always holds and 
is not particular to the steady state. The same analysis applies to the muscle isoform: 



h + f2 + h + U = fo, (7) 

where fo is the total concentration of the muscle isoform. 

At steady state the rates of change of all variables are zero, so the derivatives in 
equations iPajl -pH ]) are zero. This implies that f i = f 2 = f 3 = f 4 = —a. There are five 
unknown quantities in the model: the four reaction intermediates involving LDH-H, 61, 
62, 63, and 64, and the lactate concentration, X2- (The other metabolite concentrations 
and the metabolic flux are assumed known). The five equations required for a unique 
solution are provided by the heart isoform conservation law and the definitions of the 
reaction rates: 

ei + 62 + 63 + 64 = 60, (8a) 
kiXiCi — A;_i62 = —a, (8b) 
^22^262 - A;_2e3 = -a, (8c) 
A;363 - A;_3?/ie4 = -a, (8d) 
A;464 — /c_42/2ei = —ex. (8e) 



Equations (jHa|) . (j8b|) . (j8d|) . and (jSe|) form a linear system that determines ei, 62, 63, 
and 64. The solutions for 62 and 63 can then be substituted into (|Hcj) to determine X2. 
The result is 

ao6o - aia . 
^2 = T , (9) 
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where 



ai = A;_iA;_2^4 + k^ik^k^ + A;_2^i^42^i + kik^k^xi + A;_iA;_2fc-3?/i 

+k_2k-3kixiyi + k_ik_2k-m + k^ik_^ksy2 + k_ik^-ik^^yiy2 

+k^2k-3k-^4yiy2, 
bo = kik2k3kiXi, 

hi = k2k3k4 + ^1^2^32^! + kik2k^Xi + k_3kik2Xiyi + A;_4/c2/c32/2 
+A;_3A;_4A;2yi2/2- 

Note that when a = we recover the equihbrium relationship 

X1X2 k-ik-2k-3k-4 



yiy2 kik2ksk4 



Keg, (10) 



where i^eg is the equihbrium constant. In other words, when the metabohc flux is zero 
the system approaches an equilibrium steady state, which is independent of the total 
enzyme concentration. 

The concentration of lactate is a function of the metabolic flux, a, as shown in fig- 
ure |S] When the metabolic flux is positive (lactate is being removed from the system) 
LDH-H produces a higher steady-state lactate concentration than an equal concentra- 
tion of LDH-M. In our calculations we assumed concentrations of NADH, NAD+, and 
pyruvate equal to 0.97 /iM, 0.5 mM, and 99.4 fiM respectively. However, this quali- 
tative trend — higher steady-state lactate concentration produced by LDH-H — remains 
unchanged as long as the pyruvate concentration lies between 40 and 2 mM. 

Outside of this range the steady-state solution is infeasible (at steady state one or 
more of the metabolites has a concentration less than zero) or LDH-M produces more 
lactate than LDH-H. This interval was found by varying NADH and NAD+ indepen- 
dently over a wide range of concentrations (0.1 fiM- 10.0 mM). Previous measure- 
ments of cellular pyruvate concentration have found a range 0.08 - 0.3 mM, although 
the concentration of pyruvate has been found as high as 0.7 mM in skeletal muscle 
under tetanic conditions dStambaugh fc Post, 1966) [Boyer, 1975| Tilton et g/., 19911 



Lambeth fc Kushmerick, 2002D . Therefore, we predict that in most cells LDH-H pro- 



duces a higher steady-state lactate concentration than an equal amount of LDH-M. 

Our result implies that LDH-H produces a greater concentration of lactate in the cell 
then LDH-M does when the metabolic flux is positive. Previous work on the kinetics 
of homolactic fermentation has arrived at the opposite conclusion. The difference 
between our work and previous studies dStambaugh fc Post, 1966|) is that we have 
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focused on the steady state of the reaction. Previous analyses have adopted a Michaehs- 
Menten approach, examining the behavior of the reaction far from steady state and 
under in vitro conditions. In Michaehs-Menten theory the only possible steady state 
is the equilibrium state. Our model includes equilibrium as a special case when the 
metabolic flux is zero (see figure El). However, when the metabolic flux is nonzero the 
concentrations of the enzymes and the reaction mechanism contribute to determining 
the steady state. Thus at steady state LDH-H results in more lactate than LDH-M does 
when the metabolic flux is positive. This relationship is reversed when the metabolic 
flux is negative. 

Predicted change in lactate concentration after MAP kinase activation 

Here we predict how activating the MAP kinase pathway affects the steady-state 
lactate concentration. Treating K562 cells with PMA activates the MAP kinase path- 
way. Conversely, U0126 is a downstream inhibitor of MAP kinase signaling which acts 
on MKKl and MKK2 (Gross et ai, 20001. The array data show changes in the ex- 



pression of both LDH isoforms after activation of the MAP kinase pathway (figure QJ. 
We used the relative abundances from our array data to model three conditions: {i) 
control (cells are untreated), (ii) MAP kinase active (PMA treatment), and {Hi) MAP 
kinase partially suppressed (PMA+U0126 treatment). 

We used the full model (figure |21), with both isoforms present. The array data 
determined the relative concentration of each LDH isoform. The concentrations of the 
metabolites and enzymes in homolactic fermentation were not measured in K562 cells. 
However, metabolite and enzyme concentrations are available for muscle cells and we 
used these data as reference values in our model. K562 cells are unlikely to have an 
identical metabolic state to muscle cells, so we performed a robustness analysis to de- 
termine how our results depend on the reference concentrations used. Metabolite con- 
centrations were sampled over two orders of magnitude about the reference values. We 
verified that our results remain qualitatively unchanged over this range. We used a to- 
tal LDH concentration of 3.43 /iM QMulquiney fc Kuchel, 2003| ) in the control condition 
and assumed that the mRNA expression data directly predict protein concentrations. 
Therefore, in the control model the total concentrations of LDH-H and LDH-M were 
1.98 /jM and 1.45 /iM. For the PMA-treatment condition, the total concentrations of 
LDH-H and LDH-M were 1.83 /iM and 0.90 /iM, while for the PMA+U0126-treatment 
condition, the total concentrations of LDH-H and LDH-M were 1.27 /iM and 1.37 /iM. 
The steady-state concentration of lactate was determined numerically as described in 
Methods, with the steady-state concentrations of NAD+, NADH, and pyruvate set 
to 0.5 mM, 0.97 /iM, and 99.4 /iM, respectively, and a metabohc flux of 10.0 /iMs~^ 
( Lambeth fc Kushmerick, 2002D . 
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Our model predicts a decrease in lactate concentration for both the PMA and 
PMA+U0126 conditions (figure EI)- The expression data show that activating the 
MAP kinase pathway with PMA down- regulates LDH-H and LDH-M (figure H)). As a 
result of these changes, the model predicts a decrease in the cellular lactate level from 
37.1 ^M to 33.2 fxM. Treatment with PMA+U0126 also down-regulates both LDH 
isoforms (figure H]), leading to a predicted lactate concentration of 33.5 ijM. 

Robustness analysis 

Here we demonstrate the robustness of the model predictions. The experimentally 
observed changes in LDH-isoform mRNA predict a decrease in lactate concentration 
after PMA or PMA-I-U0126 treatment, independent of the precise values of metabo- 
lite concentrations used in the model. The values of NADH, NAD+, and pyruvate 
concentrations can vary among different cell types and growth conditions. While the 
quantitative predictions depend on the values of the metabolite concentrations, the 
qualitative result is robust. 

We used each metabolite reference value as the mean of a sampled distribution. We 
chose a lognormal distribution to describe the population of metabolite concentrations 
because samples vary by large fractional amounts ( |Boyer, 1975| ). The standard devia- 
tion is one order of magnitude: approximately 68% of the samples are within 0.1-10.0 
times the mean concentration. We note that the lognormal distribution is nonnegative, 
guaranteeing nonnegative concentrations. 

We examined 10^ randomly sampled parameter sets (figure |3)). For each set of 
parameters, we calculated the lactate concentration under three conditions: control, 
treatment with PMA, and treatment with PMA+U0126. Due to our broad sampling 
of parameter space, we found some parameter sets where the model cannot reach a 
physically valid steady state (some parameters result in a mathematical steady state 
with negative lactate concentration). Thus we excluded from our analysis parameter 
sets that resulted in an invalid steady state for any of the three conditions. This exclu- 
sion did not significantly alter the lognormal distribution of parameter values (figure 
[7j). We also checked the consistency of the simulations by comparing the distribution of 
the first 1,000 results with the distribution of the following 9,000 results. There was no 
statistically significant difference between these two distributions by the Kolmogorov- 
Smirnov test {p < 0.01) ( |Press et ai, 1992 1 . 

The median predicted lactate concentration from the randomly sampled parameters 
(figurelHI) follows the same trend observed for the reference parameter set (figurelHI). The 
predicted lactate concentration decreases by 16.9% for the PMA treament and by 14.8% 
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for the PMA+U0126 treatment, relative to control. The systematic difference between 
the three conditions can be seen by comparing sets of simulations corresponding to 
the same parameter values. The difference histograms (figure |H1) show the decrease in 
lactate concentration from the control to the treated condition. 

The predicted lactate levels for the PMA and PMA+U0126 conditions are similar. 
This result was unexpected because the LDH H:M isoform ratios were different (1.35 
for PMA and 0.85 for PMA+U0126). The key difference between the two treatments 
and the control is the decrease in the total concentration of LDH, an observation that 
led us to study the effects of isoform switching and abundance changes separately. 

Isolating the effects of isoform switching and abundance change 

We studied the effects of independently varying the isoform ratio and LDH abun- 
dance using the robustness protocol described above (figure Ej). First, we obtained 
concentrations for pyruvate, NAD+, and NADH by sampling their distributions. The 
histograms of sampled values are similar to those shown in figure |71 (data not shown). 
For each sampled parameter set we performed three simulations: (i) control, where 
the total concentration of LDH was 3.43 fiM and the isoform ratio 1.02, [ii) increased 
isoform ratio, where the isoform ratio was 1.35 and the total LDH concentration re- 
mained 3.43 fiM, and {Hi) decreased total concentration, where the isoform ratio was 
the control value (1.02) and the total concentration was 2.47 //M. These values were 
taken from the mRNA expression data for control and PMA treatment. This analysis 
allowed us to isolate the effects of changing the isoform ratio and the total concentration 
of LDH. In each simulation, the steady-state concentration of lactate was determined 
numerically as described in Methods. 

Figure ini summarizes the effect of increasing the isoform ratio or decreasing the total 
concentration of LDH. When the isoform ratio is increased, the lactate concentration 
increases, and when the total concentration is decreased, the lactate concentration 
decreases. Increasing the isoform ratio increased the lactate concentration for every 
parameter set sampled. Conversely, decreasing the total concentration consistently 
decreases the lactate concentration. The effect of decreasing the total concentration is 
greater than the effect of changing the isoform ratio, given the experimentally observed 
abundance changes in the isoforms. For these experimental conditions, we predict that 
the total abundance of LDH is more important than the isoform ratio in determining 
the lactate concentration. 



Discussion 
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Our mRNA expression data show that activating MAP kinase signahng changes the 
ratio of lactate dehydrogenase isoforms in K562 cells. We used a mathematical model 
of lactate production by LDH to calculate the changes in steady-state lactate concen- 
tration which result from changes in LDH concentration. We assume that changes 
in gene expression predict changes in enzyme concentration. We predict that for the 
experimentally observed changes in LDH, the cellular lactate concentration undergoes 
a small but significant decrease. The robustness analysis demonstrates that our predic- 
tion holds for a wide range of metabolite concentrations. Experiments which measure 
the lactate concentration in K562 cells under different conditions (control and MAP 
kinase signaling active) can directly test our prediction. 

It is often stated in the literature that the two main LDH isoforms — LDH-H and 
LDH-M — promote different directions of the homolactic fermentation reaction. The 
idea that LDH-M favors the production of lactate and LDH-H favors the produc- 
tion of pyruvate is sometimes used to explain experimental results ( Baker et al, 1997| 
Segal fc Crawford, 1994). This interpretation of the role of the isoforms was based on 



in vitro biochemistry under Michaelis-Menten conditions, which do not necessarily ap- 
ply in vivo. Homolactic fermentation is influenced by external supply of and demand 
for the reactants and products, which means the reaction is not isolated. In addition 
to assuming an isolated reaction, Michaelis-Menten theory describes the behavior of 
a reaction in its initial stages. However, metabolic reactions in the cell are typically 
close to steady state. 

In our analysis, we focus on the nonequilibrium steady state of the reaction in the 
presence of a metabolic flux, which represents the production/consumption of NADH, 
NAD-I-, lactate, and pyruvate by other sources in the cell. Under typical cellular 
conditions, LDH-H produces a higher steady-state lactate concentration than does 
LDH-M. We therefore state that LDH-H favors the production of lactate more than 
LDH-M does. However, this does not imply that LDH-M favors the production of 
pyruvate. Both isoforms favor the production of lactate when pyruvate is supplied 
to the reaction (positive metabolic flux). If the metabolic flux is negative, LDH-M 
produces more lactate than LDH-H. 

We show that examining changes in the LDH isoform ratio alone leads to incorrect 
predictions: changes in total abundance of the isoforms must also be considered. The 
changes in LDH isoform ratio we observe lead to relatively small predicted changes in 
the amount of lactate. The changes in total concentration of LDH lead to a larger 
predicted change in lactate concentration. Taking into account total concentration 
changes, as well as changes in isoform ratios, is essential for a full understanding of the 
system. 

An important problem in systems biology is the integration of information from 
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disparate sources (|Kitano, 2002D. We describe an approach to metabolic modeling that 



incorporates three important components: (i) the use of global profiling data to identify 
an interesting problem and to guide the quantitative formulation of the model; {ii) a 
kinetic model which describes the full dynamics of the system; and {Hi) a robustness 
analysis to support the conclusions. To our knowledge, no previous metabolic modeling 
work has incorporated all of these elements. The testing of our approach on this small 
problem is a pilot study for applying the method to larger systems, where the approach 
can be even more valuable. 

In recent years genome-scale profiling has become common. Significant hurdles 
remain in the interpretation and use of these data: how can information from profiling 
be integrated to advance our knowledge? In this paper, we began with an intriguing 
connection found in the experiments: activating MAP kinase signaling changed the 
expression of LDH isoforms. We predicted changes in cellular lactate metabolism based 
on the data. The careful analysis of our profiling data allowed us to isolate interesting 
and new effects. We emphasize that our model is based on experimentally observed 
changes in LDH expression, so our results are experimentally relevant. 

The advantage of using a kinetic model is that we can describe the full dynamics of 
the system, including time-dependent behavior. Valuable information about dynam- 
ics can arise from this approach; for example, the time scale required for the system 
to reach steady state can be determined. However, use of kinetic models is more 
complicated than other approaches that do not consider the full dynamics. The major 
challenge in kinetic modeling is that many parameters are unknown. Therefore, robust- 
ness analysis is essential. The robustness analysis demonstrates that our conclusions 
do not apply only for a specific parameter set, but are true in general. This component 
of our approach addresses the fundamental problem of unknown parameters in kinetic 
modeling. 

In the future we hope to apply this approach to larger systems. Testing the method 
on a smaller system (such as the one described in this paper) is an important step in the 
development of the method. However, we note that our results illustrate the power of 
carefully analyzing small systems — surprising results can be obtained through studies 
of this type. 
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Figure 1: Changes in LDH-H and LDH-M mRNA expression after treatment with PMA 
and PMA+U0126. PMA treatment activates the MAP kinase pathway, while U0126 
is a downstream MAP kinase inhibitor. The data were obtained from experiments 
performed on Affymetrix gene chips in triphcate (see text). Each data point is the 
average of three measurements and the error bars represent the maximum and minimum 
values. Treatment with PMA reduces the abundance of both LDH isoforms, and the M 
isoform shows a larger reduction. Treatment with PMA+U0126 reduces the abundance 
of both isoforms (relative to control), but the H isoform shows a larger reduction. 
The expression levels of LDH-M are significantly different among all of the treatments 
(j9 < 0.07). The expression levels of LDH-H for the PMA and PMA+U0126 treatments 
are not significantly different from one another, but they are both significantly different 
from the expression level of LDH-H in the control {p < 0.02). 
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LDH-H + NAD ^ LDH-H:NAD 
LDH-H:NAD + Lactate ^ LDH-H:NAD:Lactate 
LDH-H:NAD:Lactate ^ LDH-H:NADH + Pyruvate 
LDH-H:NADH ^ LDH-H + NADH 
LDH-M + NAD ^ LDH-M:NAD 
LDH-M:NAD + Lactate ^ LDH-M:NAD:Lactate 
LDH-M:NAD:Lactate ^ LDH-M:NADH + Pyruvate 
LDH-M:NADH ^ LDH-M + NADH 

Figure 2: The elementary reactions of homolactic fermentation. Homolactic fermenta- 
tion is a compulsory-order, ternary reaction. All elementary reactions are reversible. 
The arrow in each elementary reaction indicates the direction of a positive reaction 
rate. 
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Figure 3: Rate constants in the kinetic model ( Borgmann et aL, 1975D 
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Figure 4: Schematic of the computational protocoL The concentrations of pyruvate, 
NAD-I-, and NADH are determined in one of two ways: either they are specified (basic 
protocol) or they are randomly sampled (randomized protocol). The total concentra- 
tions of LDH-H and LDH-M arc determined from the experimental conditions with 
both isoforms present. We determine the steady-state concentration of lactate using a 
hybrid approach. The model will relax to steady state if integrated sufficiently long, 
but Newton's method accelerates the convergence to steady state. 
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Figure 5: The predicted concentration of lactate as a function of the metabohc flux a. 
In the model, a non-equilibrium steady state is possible when pyruvate and NADH are 
added to the system at a constant rate a, and lactate and NAD+ are removed at the 
same constant rate. These terms in the model represent the production of pyruvate 
(and consumption of lactate) in other chemical reactions or transport into/out of the 
cell. As the metabolic flux a is varied, the steady-state lactate concentration predicted 
by the model changes. Note that a > means that lactate is being removed from 
the system. In each of the curves shown, only one of the LDH isoforms is present. 
The concentrations of the metabolites and the total amount of each enzyme are held 
constant (pyruvate, 99.4 /zM; NAD+, 0.5 mM; NADH, 0.97 ^M; LDH, 3.43 ^M). The 
sign of the metabolic flux determines which of the two isoforms favors the production 
of lactate. This conclusion differs from previous analyses of this reaction, which were 
performed under Michaelis-Menten conditions (see text). The experimental verification 
of this prediction is future work. 
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Figure 6: Predicted concentration of lactate for three conditions which simulate control, 
treatment with PMA, and treatment with PMA+U0126. The three conditions differ in 
the activity of the MAP kinase pathway. PMA treatment (which activates MAP kinase 
signaling) results in a small but significant decrease in predicted lactate concentration 
of 10.5%. Treatment with PMA+U0126 slightly increases the lactate level relative 
to the PMA treatment. The similar lactate levels for PMA and PMA+U0126 are 
surprising, because the LDH H:M ratio is 1.35 for the PMA treatment and 0.85 for the 
PMA+U0126 treatment. The concentrations of pyruvate, NAD+, and NADH are 99.4 
fiU, 0.5 mM, and 0.97 /xM. The metabohc flux a is 10.0 /iM s'^. 
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Figure 7: Sampled metabolite concentration values. The distribution of each con- 
centration is lognormal (see text). The mean of each distribution coincides with the 
reference concentration of that metabolite, and the standard deviation is one order of 
magnitude. 
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Treatment Change in lactate (M) Change in lactate (M) 

Figure 8: Treating K562 cells with PMA or PMA+U0126 is predicted to result in a sim- 
ilar concentration of lactate. (A) The median concentration of lactate under the three 
conditions examined in the robustness analysis (control, PMA, and PMA+U0126). 
This result is qualitatively similar to the result from the simulations using only the 
reference concentrations of pyruvate, NAD+, and NADH (figure Treatment with 
PMA or PMA+U0126 decreases the concentration of lactate. (B) Histogram of the 
predicted lactate concentration for the control simulation. This panel illustrates the 
distribution of the lactate concentration found in the robustness analysis. The median 
of this histogram is the value of the control treatment shown in A. (C) Histogram of the 
difference between the lactate concentration predicted for the control and PMA treat- 
ment conditions. All of the differences are positive, which means that PMA treatment is 
predicted to decrease the steady-state lactate concentration. After PMA treatment, the 
LDH isoform ratio increases from 1.02 to 1.35. (D) Histogram of the difference between 
the lactate concentration predicted for the control and PMA+U0126 treatment con- 
ditions. All of the differences are positive, which means that PMA+U0126 treatment 
is predicted to decrease the steady-state lactate concentration. After PMA+U0126 
treatment, the LDH isoform ratio decreases from 1.02 to 0.85. 
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Figure 9: The contrasting effects of increasing tlie isoform ratio and decreasing tfie total 
concentration of LDH. (A) Tlie median concentration of lactate under the three con- 
ditions examined (control, isoform-ratio increase, total LDH concentration decrease). 
When the isoform ratio is increased and the total concentration of LDH is held con- 
stant, the amount of lactate produced increases by a small amount. This trend is 
consistent with the behavior of the individual isoforms: LDH-H produces more lac- 
tate than LDH-M for a > 0. When the isoform ratio is held constant and the total 
concentration of LDH is reduced, the amount of lactate produced decreases. (B) His- 
togram of the predicted lactate concentration for the control simulation. This panel 
illustrates the distribution of the lactate concentration found in the robustness analy- 
sis. The median of this histogram is the value of the control treatment shown in A. (C) 
Histogram of the difference between the lactate concentration predicted for the isoform- 
ratio increase and control conditions. Every set of metabolites resulted in an increase 
in the predicted lactate concentration. (D) Histogram of the difference between the 
lactate concentration predicted for the total LDH concentration decrease and control 
conditions. Every set of metabolites resulted in a decrease in the predicted lactate 
concentration. Note that the magnitude of the change in the lactate concentration is 
approximately 10 times larger than the effect of changing the isoform ratio (shown in 
C). The influence of the total concentration of LDH on the production of lactate is 
much greater than the influence of the isoform ratio. 
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